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BioScape is a concurrent language motivated by the biological landscapes found at the interface of 
biology and biomaterials (5). It has been motivated by the need to model antibacterial surfaces, 
biofilm formation, and the effect of DNAse in treating and preventing biofilm infections. As its 
predecessor, SPiM lfl2l . BioScape has a sequential semantics based on Gillespie's algorithm Q, and 
its implementation does not scale beyond 1000 agents. However, in order to model larger and more 
realistic systems, a semantics that may take advantage of the new multi-core and GPU architectures 
is needed. This motivates the introduction of parallel semantics, which is the contribution of this 
paper: Parallel BioScape, an extension with fully parallel semantics. 

Process algebras have been successfully used in the modeling of biological systems, see lfl4l l4l[TTl. 
where they are particularly attractive, because of their ability to accommodate new objects and new be- 
havioral attributes as the complex biological system becomes better understood. However, most of the 
modeling languages lack adequate support for the design of systems in which to study complex inter- 
actions involving both spatial properties, movements in three-dimensional space, and stochastic inter- 
actions. Recently, new spatial modeling languages allowing explicit description of temporal spatial dy- 
namics of biochemical processes have been proposed (SpacePi [8 ], DCA lfT7l . LYl |[T6l . Stochsim iflOlD . 
Other agent-based platforms (9] include C-Immsim |[T5l[TTI and PathSim visualizer 031 ■ However, few 
of them support individual based, continuous motion, and continuous space stochastic simulation O, 
which are important features for modeling temporal spatial dynamics of biochemical processes accu- 
rately. To address this problem in previous work we introduced BioScape @, a language incorporating 
both stochasticity and 3D spatial attributes. 

Gillespie's algorithm produces two outputs in each iteration: 1) the next reaction R to be executed 
and 2) a slice of time t to advance the simulation clock. Since many reactions, including many instances 
of the same reaction, may be available, the slice of time t does not correspond to the time that R would 
take, but an amount of time proportional to the time it would take to execute all available reactions. In 
contrast, the parallel semantics will execute all available reactions, not just one instance of one reaction 
R, and the first challenge is then how to calculate simulated time. Reaction times can vary substantially, 
for example, some prokaryotic cell mitosis takes ten minutes, some plant cell mitosis takes about half an 
hour, while some animal cell mitosis takes about three hours. If we trigger all reactions together, how do 
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P,Q ::= | X(u) \ P \ Q | (va@r,rad).P D ::= | D,X(x) = M^ a '° FV(M) Cx 
M::= n.P[+M] u,v::= a \ b \ ■■■ \ x | y | ••■ 

K ::= delay@r | \u(v) \ lu{x) | mov E ::= | £,a@r,rad 

Figure 1: Syntax 

we advance the simulation clock? The solution we propose here consists of annotating each product of a 
reaction with a timer indicating how long that reaction will take. 

For example, if Cell — >3o Cell \ Cell means that a Cell takes 30 minutes to split, through mitosis, 
into two daughter cells, then we will annotate the two daughter cells as {{Ce//}} 30 and {{Ce//}} 30 . As 
time lapses, the timer will be reduced, and when reaching {{Ce//}}°, both cells will be available for new 
reactions. 

In Fig. Q] we define the syntax of the calculus, which slightly simplifies the syntax of 121 in order to 
avoid decorating semantic processes with shapes, as defined at page [103] 

We assume a set of channel names, denoted by a, b, and a set of variables, denoted by x, y, with 
subscripts or superscripts, if needed. As usual, aisa\,... ,a n , and similar for x. The empty process is 
0. By Xiu) we denote an instance of the entity denned by X. The actual parameters of the instance may 
be either channel names or variables, in case the instance occurs in a definition. The process P \ Q is the 
parallel composition of processes P and Q. By (va , @r,rad).P we define the channel name a with two 
parameters r and rade M>o within process P; the parameter r is the stochastic rate for communications 
through channel a and rad is the communication radius. The radius is the maximum distance between 
processes in order to communicate through channel a, and the reaction rate determines how long it takes 
for two processes to react given that they are close enough to communicate. 

The heterogeneous choice is denoted by M, where K.P [+ M] means K.P \ K.P + M. Choices may 
have reaction branches and movement branches. The reaction branches are probabilistic (stochastic), 
since reactions are subject to kinetic reaction rates, while the movement branches are non-deterministic, 
since the movement of instances of entities is always enabled, provided there is enough space. The 
prefix 71 denotes the action that the process n.P can perform. The prefix delay@r is a spontaneous and 
unilateral reaction of a single process, where r is the stochastic rate. The prefix lu denotes output, and 
the prefix ?w denotes input. The prefix mov moves processes in space according to their diffusion rate 
(ft)) (see below). We use standard syntactic abbreviations such as % for 7T.0. 

We denote by D a global list of definitions. The equality X(x) = M^- W,<J defines entity X with formal 
parameters x, to be the choice M with geometry ^ , CO, G, specifying a movement space t, , a step ft), and 
a shape a. The choice M describes the behavior of X with a choice of prefixed processes. The selection 
of one of the choices depends not only on the available interactions with other processes, but also on 
the available space. The movement space § is a set of point coordinates in the global coordinate system 
defining a volume. Intuitively, X can move within <^ . The step CO € M>o, is the distance that X can stir 
in a movement, and it corresponds to the diffusion rate of X; a is the three-dimensional shape (sphere, 
cube, etc.) of X. The movement space for the empty process is everywhere, the global space, and 
its movement step is by default. The entity variable X can be defined at most once in D, and the free 
variables of P, denoted by FV(P), must be a subset of the variables x. We also write X(x) = (k.7i' .Pp' a ' a 
as short forX(x) = (tc.Y (x))^ m < a and Y(x) = (n' .P)^' a ' a . 

We use E to range over environments of channel name declarations. By a@r, rad we declare channel 
name a with reaction rate r and reaction radius rad. A channel name a appears at most once in E. 

Consider the following simple example of a bacterium Bac, that can either move or divide into two 
daughter cells. Bac is defined with movement space movB, movement step stepB, and shape shapeB. 
Intuitively, bacteria can move within movB, with non-deterministic steps of length stepB, and the shape 
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S.Loc _ . S.Loc.Par 

P=Q SXOC.INU Ml (shape(P)) UM2(shape((2)) = /I (shape (P | Q)) 

= {2}^ (va@r,rad).{P} M = {(Vfl@r,rad).F} M Mm, I {G}/i 2 = i P I 2>M 

S.Nu.Par 



S. Nu.Com 



' f n(B) 



(va@r,rad).(W>@r',rad').A = (vfe@r',rad').(va@r,rad).A ((va@r,rad) A) | fl = (va@r, rad).(A | B) 

Figure 2: Structural Equivalence of Spatial Configurations 

shapeB is at all times contained within movB. The prefix mov represents a non-deterministic move- 
ment of length stepB, whereas delay@1.0. (Bac() I BacQ) represents mitosis, the division of a 
bacterium into two daughter cells: Bac() I Bac (), and the del ay@l . prefix is used to model the fact 
that division is not an instantaneous reaction. 

Bac() = (mov.Bac() + delay® 1.0. (Bac()|Bac())) movB - ste P B - sha P eB 
A run-time system is represented by a parallel composition of entity instances (without free variables) 
each with its shape, and located in some positions of a global frame. We define the shape of processes 
inductively as follows: 
shape(O) = shape(X(a)) = a if X(x) = M$> m > a € D 

shape (( va@r,r ad) .P) = shape(P) shape(P | Q) = shape(P) W shape(2) 

where gives a shape obtained by composing two shapes trough juxtaposition. For different applications 
we can choose suitable functions to realise y, we only require y to be a commutative and associative 
operator, i.e. 0^02 = 0^ &i and (<Ti y 02) y a 3 = <Ti y (<r 2 y (73). 

We use n to denote a map which applied to a shape locates it in the global space, by putting its 
barycentre at a fixed point, orienting the shape, and possibly modifying it. So ju(shape(P)) computes 
the space occupied by a process P in the global coordinate system. Processes may also share channels 
for communication. Spatial configurations, denoted by A, B, ... are defined as follows: 

A,B ::= {P}^ | A\B \ (va@r,rad).A 
Structural equivalence on configurations is defined in Fig. [2j omitting the rules for associativity and 
commutativity of | and +. Parallel composition has neutral element {0}^ for any p. Rule S.LOC uses 
the standard structural equivalence of Pi-calculus processes. The premise of rule S.LOC.PAR assures 
that the two equivalent processes occupy exactly the same space. In rule S.Nu.Par, f n is a function that 
returns the set of free channel names of a configuration. 

The (parallel) operational semantics of BioScape is based on two auxiliary reduction relations: a 
stochastic relation, E h A^B, for reactions such as synchronisation and delay, defined in Fig. [3j and a 
non-deterministic (non-stochastic) relation, A— >B, for geometric transformations, in our case movement, 
defined in Fig.H] Notice that reduction axioms (SR.DEUAY, SR.Com, NR. Move) only involve entities 
(X(a)), and entities evolve according to one of the choices in their definitions. In rules SR.DEUAY, 
SR.COM and NR. MOVE, there is no check of whether the entities of the resulting process have enough 
space, since this check is done in the parallel reductionrules PR.STOC, and PR. MOVE of Fig. [5] In 



SR.Delay SR.Str 

X{x) = (delay@-r.P[+M])^ ma ieD A=A f EhA'^B 1 B 1 = B 

E h {X (a) } M ^{P[a/x] } M E h A^B 

SR.COM 

X{x) = {\a(b).P [+ M])^' e D Y{y) = {7a(z)-Q [+ N}f M '- & € D i±a(ji,fi') < rad 
£, fl @r,radh {X(c)} M | {Y(d)}^^{P[c/x}}^ | {Q[d /y][b / z]} ^ 
Figure 3: Stochastic Reduction Relation 
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NR. Move NR.Str 

fl' = translate fi'(a) C £, X(x) = (mov.P [+ M])^' a> ' a e D A=A' A'^B' B' = B 

{X(a)}^{P[a/x]}^, A^B 
Figure 4: Non-stochastic Reduction Relation 

particular, a stochastic (non-stochastic) redex is stuck, if there is not enough space for its reduct in the 
configuration. Therefore, the evolution of systems in parallel BioScape produces configurations in which 
space is consistent. 

Fig. [3] defines the stochastic reduction relation of BioScape, E h A-^-B, where r is the rate of the 
channel used for synchronization or delay. We write dis(jU,jU') for the distance between the origin of 
p. and the origin of p' . In rule SR.COM the condition dis(p,p') < rad ensures that located processes 
{X(c)} ll and {Y(d)} ll i are close enough to communicate through channel a. The non-stochastic re- 
duction relation of BioScape, A— >B, is defined in Fig. @] By translate(o),/i) we denote the function 
that randomly generates a new map p', using the movement step CO and the old map p. The condition 
ju'(ct) C <§ of rule NR. MOVE ensures the new located process {P[a/x]}^i is within the movement space 
t, of X (see previous remark about not checking if the entity moves to an empty space). 

For stochastic reductions we compute the duration of the reduction, based on the exponential distri- 
bution associated with the propensity of the reduction. Since reductions may have different durations, 
we introduce timed configurations, {{A}}", meaning that, after a time n, this configuration will be A. 
We extend structural equivalence to timed configurations by adding that {{A}} = A, and A = B implies 
{{A}}" = {{#}}". With the metavariables F, and G we denote either spatial configurations or timed con- 
figurations (extended configurations), i.e., 

F,G::=A | {{A}} n \ F\G | (va@r,rad).F (n > 0) 

We define a reduction strategy that given the whole configuration, first moves all the processes that 
can be moved, and then executes all the stochastic reductions that can be executed, omitting only reduc- 
tions which would lead to overlaps, i.e. configurations where some entities occupy the same space. Both 
non-stochastic and stochastic reductions are applied in parallel. For this purpose, we define multi-hole 
contexts C by the following grammar: 

C::=F | [] | C\C | (vx@r,rad).C 
Congruence on multi-hole contexts is naturally induced by the congruence on configuration, associativity 
and commutativity of the parallel operator, and standard rules for v restrictions similar to S. Nu.COM 
and S.Nu.Par. Given this congruence any multi-hole context, C, may be written in a canonical form. 

That is, there is C', C = C' such that C' = V\ V„.Fi | • • • | F m | | • • • | [ ], where V,-, 1 < i < n, is an 

abbreviation for va^r^rad;, and for all j, 1 < j <m, Fj = {{A}}" for some A, and n, or Fj = {P}^ 
for some P, and p. We say that a\ @ri,radi, . . . ,a„@r n ,rad„ is restr(C). In the following we assume 
that multi-hole contexts are always in canonical form. 

As already mentioned, our reduction strategy avoids spatial overlaps. In particular for moving re- 
ductions we have to ensure that moves and reshaping are compatible with the available space, that is 
after moving no entity overlaps with another entity. For stochastic reductions we have to assure that the 
created entities have their space. To this aim we define the space of a configuration, and a predicate that 
says whether a configuration does not have any overlapping entities. 

Let space(F) be a function on configuration F that returns the space occupied by its processes 
located in the global frame defined as follows. 

space({P} p ) = p (shape (P)) space({{A}}") = space(A) 

space(F | G) = space(F) U space(G) space((vfl , @r,rad).F) = space(F) 

We say that a configuration F is OK if the various entities in F do not overlap, that is: 
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PR.Move PR. Timed 

Fi->Gi (l<i<p) C[Gi] ■ • ■ [G p ] OK mv n = min{n t | 1 <i<p] C is untimed 

C[F l ]---[F p ]^C[G l }---[G p ] C[{{A! }}»']■ • ■ [{{A p }}"»] - C^A, }}'" • ■ ■ [{{A, }}»<-»] 

PR.Stoc PR.Conf 
restr(C) hA;-^; n; = T(r/,C/[A,]) (1 < i < /») C[Ci(BiJ] ••• [C p [Bp]] OK w F^F l -*F 1 ^F l 

C[C, [Ai]] ■ • ■ [C p [A p \] — C[d [{{Si}}">]] ■ • ■ M^Fii F^F' 

Figure 5 : Parallel Reduction Relation 

{P} M OK A OK => {{A}}' 1 OK FOK ^ (vi@r,rad).F OK 

F OK A G OK A space(F) n space(G) =0 F \ G OK 
With the notion of OK configuration we define two notions of well-formedness of configurations. The 
first notion is to be used for parallel move reductions and the second for parallel stochastic reductions. 
Theses notions are to be used to enforce (z) the fact that only reductions that have enough space for 
their reduct are allowed, and (//) that we want maximal parallelism, that is any "extra" movement or 
transformation would produce an overlap. In order to formalise this we first need to single out the sets 
9? mv and 9?^ of movement and stochastic redexes, i.e. we define: 

• % nv = {{X(a)}^ | X(x) = (mov.P+M)^ w ' a GO}, 

• % st = {{X(a)}^ | X(x) = (delay@T.P + M)^' (0 ' a G D}U 

{{X(c)} M | {Y(d)}^ | X(x) = (\a(b).P + M)^ a ' a eD&Y(y) = (7a(z).Q+N)^'' (0 '' a ' e D 
& dis(jU,jU') < rad} where a@r,rad is the declaration of channel a. 

We extend the syntax of configurations by allowing underlined extended configurations, defined by: an 
underlined extended configurations is a configuration in which some spatial sub-configurations may be 
underlined. Underlined configurations are the tool we use to define maximal parallelism. We can then 
define: 

Definition 1. ( i) An extended configuration F is OK mv if F is OK and F = C[A] with A not underlined 
and A 6 9t mv and A^B imply C[B] not OK. 

(ii) An extended configuration F is OK if ifF is OK and F = C[A] with A not underlined and A € 3i st 
andA^B imply C[B] not OK. 

As a last notion, we say that a context C is untimed if it does not contain timed configurations. 

We are now able to explain our parallel reduction strategy, whose rules are given in Fig. [5] The first 
three rules deal respectively with parallel movements, timed reductions, and stochastic reductions, while 
the fourth rule maps extended configurations into extended configurations by applying first the parallel 
movements, then the stochastic interactions, and finally by advancing the time of the minimum required 
to complete one or more interactions. In this way at the next iteration there would be new entities to be 
moved and/or stochastically reduced. 

The condition of obtaining an OK, m , extended configuration in rule PR. MOVE assures that all pos- 
sible moves in C[F\] ■ ■ ■ [F p ] which do not cause overlaps have been done in the reduction. Similar effect 
is produced by the conditions that the extended configuration is OK sf and that the context is timed in 
the following two rules, respectively. Rule PR.STOC prescribes that the time of a stochastic reaction 
depends (through the function t) on the rate of the reduction and on the number of available reactants. 
The outer context C is a multi-hole context, while the context C, of the redex Aj is a single hole context 
capturing the surrounding environment that influences the speed of the reduction. We could incorporate 
a counting function keeping track of the available reactants in the communication range (in a way similar 
to what is done, e.g., in El|6l). 

Examples, results of simulations, comparisons with related papers and discussions can be found in the 
full version of this papers available at http : //www.di .unito . it/~dezani/papers/cdgsst .pdf 
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